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We present a Direct Statistical Simulation (DSS) of jet formation on a /3-plane, solving for the 
statistics of a fluid flow via an expansion in cumulants. Here we compare an expansion truncated at 
second order (CE2) to statistics accumulated by direct numerical simulations (DNS). We show that, 
for jets near equilibrium, CE2 is capable of reproducing the jet structure (although some differences 
remain in the second cumulant). However as the degree of departure from equilibrium is increased 
(as measured by the zonostrophy parameter) the jets meander more and CE2 becomes less accurate. 
We discuss a possible remedy by inclusion of higher cumulants. 



PACS numbers: 47.27.wg, 47.27.eb, 92.60.Bh, 92.10.A- 

Jets are relatively narrow bands of fast-flowing fluid 
that move coherently in one direction. They are ubiqui- 
tous in nature, found in Earth's oceans and atmosphere, 
the outer layers of gas giant planets, the interior of stars, 
and laboratory experiments with fluids and plasmas pQ. 
Jets play an important role in the dynamics of these 
fluids and it is therefore important to understand the 
mechanism(s) that govern their formation. It is often the 
case that jets are driven by energy input at small spatial 
scales; the question then is how this energy is transferred 
into large scale coherent motion. Two competing mech- 
anisms have been proposed, both of which rely on the 
interaction of turbulence and rotation. The first involves 
the scale-by-scale transfer of energy known as the inverse 
cascade [2]. Large-scale vortices are known to be gen- 
erated by this mechanism. The other mechanism relies 
upon direct transfer of energy to the largest scales. It 
is hard to disentangle these two mechanisms in experi- 
ments and in simulations. Direct calculation of statistics 
and quasilinear direct numerical simulation (DNS) cal- 
culations have demonstrated that jets can be formed by 
the direct mechanism, not relying on an inverse cascade 
|3HS1- Jets control mixing, and in the case of the Earth, 
the position of the storm tracks; therefore improved un- 
derstanding of jets will lead to better models of climate 
and other complex systems. 

Non-equilibrium statistical mechanics can be used to 
focus attention on universal aspects of fluids in motion. 
Isotropic, homogeneous, turbulence is at present beyond 
the reach of a complete statistical theory. By contrast, 
inhomogeneous flows (and their magnetohydrodynamic 
counterparts) such as jets may be accessible to Direct 
Statistical Simulation (DSS), that is, methods that solve 
directly for the statistics of the flow [7]. DSS offers the 
possibility of a deeper understanding of mechanisms of 
fluid motion, as well as a practical speed-up in obtaining 
statistics [5]. In the limit of small driving and dissipa- 
tion, equilibrium statistical mechanics is a powerful tool 



for understanding quasi two-dimensional flows (for a re- 
view see Ref. [8 ). Away from equilibrium, Stochastic 
Structural Stability Theory (SSST) [I [U [10] is one ap- 
proach that has been explored to understand the forma- 
tion and maintenance of jets. Here we instead investigate 
systematic expansion in equal-time, but spatially nonlo- 
cal, cumulants of the flow. When truncated at second 
order, the cumulant expansion (denoted CE2) is closely 
related to SSST [6] but it is only the starting point for a 
perturbative expansion in higher cumulants. 

This paper examines the accuracy of DSS at the CE2 
level of approximation as a representation of the statistics 
of turbulent flows driven far from equilibrium. CE2 in- 
cludes the interaction of mean flows with eddies to drive 
eddies and the interaction of eddies with eddies to drive 
mean flows, but removes the interaction of eddies with 
eddies in the evolution equation for the eddies [11] ; an in- 
teraction that has been termed the "EENL" (eddy-eddy 
nonlinearity) by Srinivasan and Young [6 . It has been 
argued [HI [12] that CE2 is an exact representation in the 
quasi-equilibrium limit, but the domain of validity of such 
a truncation as the system is driven further from equilib- 
rium remains largely untested. We conduct numerical ex- 
periments to investigate the accuracy of DSS at CE2 for 
systems removed from quasi-equilibrium by considering a 
model problem of the driving of jets by small-scale forcing 
on a /3-plane. This system has been studied extensively 
within the framework of DNS in both the fully nonlinear 
and quasi-linear regimes (see for instance Refs. [6l H3])- 
Although this model is the simplest that includes all the 
requisite features for our purposes, i.e. anisotropy, non- 
trivial long-range correlations and mean flows, we note 
that it is a particularly rigorous test of these statistical 
methods in that it is stochastically (rather than deter- 
ministically) driven and translationally invariant in two 
directions, with only the emergence of jets spontaneously 
breaking the latitudinal symmetry [6] and leading to in- 
homogeneity. We return to this in the discussion at the 
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end of the paper. 

The /3-plane we use is periodic in both the x (longi- 
tude) and y (latitude) directions, with the domain of 
size 2tt x 2ir. The motion of the incompressible fluid 
is damped by a single friction n and also by small-scale 
dissipation that absorbs structures at the finest scales. 
(Some models examined in Ref. [10] have friction damp- 
ing the fluctuations that is 10 times greater than that 
slowing the zonal mean flow.) The fluid is driven by 
random (stochastic) forcing rj. This type of stochastic 
forcing is widely used as a simple model of small-scale 
processes that inject energy into the fluid, with the small 
and fast scale acting as a random influence on the large 
and slow scales [T4HT6] . The time-evolution of the rela- 
tive vorticity £ = z • (V x u) is given by for example Ref. 

d t c + J{i>, c) + /3d x ^ = -< + ^v 2 c + v(t), (i) 

C = V 2 </>, (2) 

where J(a, b) = d x a d y b — d y a d x b. Here if; is the stream- 
function and the fluid velocity u = (u, v) = {—d y i/j, d x ip); 
we have set the deformation radius of the flow to be in- 
finite for simplicity. The forcing is chosen to be random 
with a short (but non-zero) renewal time (0.1 < r t < 1) 
and to be concentrated in a spectral band of wavenum- 
bers between k m i n < \k x \,\k y \ < k max (and for these runs 
kmin = 7, k max = 10). The amplitude of the forcing is 
chosen from a Gaussian distribution with standard devi- 
ation CLrj. This is a popular choice of forcing; a detailed 
discussion of the role of forcing in DNS of such problems 
is given in Ref. [T3] . 

The dynamics of this system has been extensively stud- 
ied since the pioneering calculations of Rhines [18 . He 
demonstrated how correlations between nonlinear Rossby 
waves could lead to the generation of zonal flows and 
identified the scale at which zonal flows become impor- 
tant in mediating the dynamics of the propagating non- 
linear Rossby waves (see e.g. Ref. [19 ]). This scale, now 
termed the Rhines scale, is given by Lr = ([///?) 2 , where 
U is the rms velocity of the flow, and occurs when the 
second and third terms of equation ([!]) are comparable. 
There has been much research into the importance of this 
lengthscale for the ultimate latitudinal scale of jets (see 
e.g. Ref. [20 and the references therein), but it is also 
becoming clear that the dynamics of zonation is also con- 
trolled by another length scale L e , originally introduced 
in Ref. [21], which measures the intensity of the forcing 
relative to the background potential vorticity gradient. 
For the simple /3-plane model L e = (e//3 3 )s where e is 
the energy input rate of the stochastic forcing rj{t). 

The ratio of these two length scales is now believed to 
play a critical role in determining the strength and sta- 
bility of jets [T7J|22], for cases where the same damping 
is applied to the mean flows and the turbulent fluctua- 
tions. This local measure, termed the zonostrophy index, 




FIG. 1. Snapshot density plot of vorticity together with 
zonal mean vorticity profile of jets found by DNS. The param- 
eters are k = 10 -3 , v — 10 -4 , ft — 16. For these parameters 
Rp = 3. 




FIG. 2. Hovmoller diagrams of zonal mean relative vorticity 
versus time from DNS simulations, (a) Parameters as for 
Fig. [I] with Rp — 3; (b) Parameters as for (a) but ft = 8 and 
Rp — 2.8; (c) Parameters as for (a) but k — 10 -2 and Rp = 
1.75. A running time- average commences at the midpoint of 
each diagram. 

is given by Rp = L R /L e = U 1 / 2 /3 1 / 10 /e 1 / 5 . In general, 
if the zonostrophy index is large {Rp > 6) then strong 
stable jets are found, whilst for small Rp the jets are 
weaker, meander more and no staircase is formed [T3] . 
The zonostrophy index is therefore a measure of how 
far the system is driven out of equilibrium. We note 
here that Rp can also be written in terms of the ratio 
of an advective time on the Rhines scale to a dissipa- 
tive timescale {Fp = kLr/U) i.e. Rp = Fp 1 ^ 5 . Hence 
the quasi-equilibrium limit given by Rp — >• 00. Recent 
estimates have put Rp > 10 for flows on the surface of 
Jupiter, whilst Rp ^ 2 for oceanic jets [T7] . 

DNS is performed using a pseudo-spectral scheme op- 
timised for use on massively parallel machines [23 . For 
these simulations we utilise resolutions of up to 2048 2 . 



The forcing is applied at moderate scale (with r t a^ = 
0.01 for all calculations) and the system is evolved from 
a state of rest until a statistically steady state is reached. 
Fig. [T] shows a snapshot of the vorticity and zonally av- 
eraged vorticity for a state with 3 zonal jets coexisting 
with a small-scale turbulent flow. For this calculation 
Rp = 3.0 and so the jet is well removed from the quasi- 
equilibrium limit. We note that this limit is exceed- 
ingly difficult to simulate in DNS, requiring extremely 
long integrations. Nonetheless the Hovmoller diagram in 
Fig. [2|a) of the (£, y) dependence of the mean flow to- 
gether with a running time average calculated from the 
midpoint of the calculation shows that the zonal flows do 
not meander too much in space and well-defined averages 
can be obtained — though we note that the dynamics 
needs to be integrated for a long time to obtain mean- 
ingful flow statistics. Fig. |5|b) shows the correspond- 
ing diagram when Rp has been reduced to 2.8, which is 
achieved here by lowering (3. For this case, even further 
from quasi-equilibrium, the jets are still relatively steady, 
but the Rhines scale has changed sufficiently that now 
only two jets fit in the domain. For the final case of 
Fig. [2|c) Rp has been reduced to Rp = 1.75, achieved 
by increasing the friction. Here the jets meander signifi- 
cantly and merge showing a large degree of spatiotempo- 
ral variation. In this respect they have the characteristics 
of observations and simulations of oceanic jets [24j [25] . 
At different times there appears either four or five jets, 
but on average there are five jets. Because of the tem- 
poral variability, the time average of the jet velocity is 
much smaller than the instantaneous jet speeds. 

Expansion of equal-time cumulants at order CE2 is 
straightforward for Eq. [I] Let r = (x,y) and r' = (V, y') 
and then adopt a Reynolds decomposition for the vari- 
ables by setting £(r) = (C( r )) + C'( r )? where the angle 
brackets imply either an ensemble average or an aver- 
age over longitude (x) that are equivalent for this for- 
malism. We define the first cumulant for the vorticity as 
C C = (C( r )) = c c(y)- It i s convenient to define the first cu- 
mulant for the streamfunction as c^(r) = (*p(r)) = c^(y) 
where the relationship between the first cumulants is 
given by cq = dyyC^. We may then define the second 
cumulants as follows: c ( ^(r,r / ) = (C( r )C( r ')} an d note 
that for this system depends on the two local latitudes 
and the difference between the longitudes £ = x — x' , i.e. 
c^(r,r') = C£(;(y,y',£) [TT]. Corresponding definitions 
arise for the derived second cumulants and c^, i.e. 
c^(r,r') = (^'(r)C'(r')) = tyc(y,y f ,£) and similarly for 
c Cip(y->y' ->C)' With these definitions, it is straightforward 
to derive the equations for cumulant hierarchy, truncated 
at second order, as 

dtc c = -(d y + d y >) d^^l, - kc c + vd% y c c . (3) 
d t c a d y c^d^c a - d y (c c (y) - f3y)d^ c 
- dy'C^c a + d y >(c c (y') - (3y f )d^ 
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v (V 2 + V' 2 ) c cc - 2kc cc + r. 



(4) 



FIG. 3. Hovmoller diagrams of relative vorticity from 
CE2. (a) Parameters as for Fig. |5Ja). (b) Parameters as 
for Fig.[2jc). 



Here T is the covariance matrix of the stochastic forc- 
ing that enters into the equation for the second cumu- 
lant as a deterministic source term localised at the same 
wavenumbers as for the DNS [5] and with an amplitude 
ar that is given ar = r t a^. Eqs. [5] and 4] constitute 
a realizable closure and in the absence of damping and 
forcing, conserve linear momentum, energy, and enstro- 
phy. Eqs. [3] and [4] are integrated forward in time using a 
pseudo-spectral scheme that utilises a mixed integrating 
factor/ Adams Bashforth timestepping scheme; details of 
the method are included in a subsequent paper. The in- 
tegrations were performed at a typical spectral resolution 
of 16 x 128. Restricting \k x \ < 16 does not amount to 
a further approximation beyond CE2, because, for this 
problem, at level CE2 only modes with zonal wavenum- 
bers less than those of the stochastic forcing are excited. 

Recall that for Rp large the system is in quasi- 
equilibrium, dominated by strong jets, and CE2 should 
provide an accurate representation of the statistics of the 
fluid flows. A typical evolution of the cumulant system is 
shown in a Hovmoller diagram in Fig.ga). After some 
initial transients where jets are driven with a relatively 
small latitudinal extent, broader jets emerge via a series 
of jet mergings. Similar jet-merging behaviour has been 
observed both in DNS and in the strong jet simulations 
of SSST [9]. The system eventually reaches a statistically 
steady state, represented by a simple fixed point of the 
cumulant equations. The calculations were repeated at a 
range of Rp and compared with the (zonal and time aver- 
ages of the) DNS solutions described earlier. Fig.[4]shows 
comparisons of the zonal velocity in the jet from DNS av- 
eraged over both x and time with that achieved from DSS 
at CE2 for Rp = 3.0 and Rp = 2.8. It seems as though 
the agreement in the first cumulant at these levels of lack 
of equilibrium is good; CE2 reproduces both the correct 
number of jets and their strength; although CE2 slightly 
overestimates the average jet strength — a characteristic 
in common with quasi-linear DNS of jets [6]. However 
close examination of the second cumulant reveals that 
CE2 struggles to reproduce the cross-correlation patterns 
(or teleconnections) from DNS for these parameters. The 
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FIG. 4. Comparison of mean zonal velocity from DNS 
(dashed lines) and CE2 (solid lines) for parameters as in 
Figs.gja), (b), and (c). 
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FIG. 5. Second cumulant c^^ as calculated from DNS 
(left) with R/3 = 3.0 and the corresponding CE2 solution 
(right). Cross correlation with respect to a test point located 
at (tt,4.7). 



left panel of Fig. [5] shows the second cumulant as accumu- 
lated from the DNS solution of Fig. [TJ The figure shows 
the cross-correlation of the vorticity statistics with re- 
spect to a test point located at a longitude of tt and a 
latitude at the maximum of the large-scale vorticity (at 
the wings of the jet). The second cumulant is localised 
in latitude, with some recurrent correlations occurring 
on the jet spacing, whilst the structure in longitude has 
both wavenumber k x — 1 contributions and contribu- 
tions from the scale of the spectral forcing. Examination 
of the spatio-temporal dynamics of the system indicates 
that the k x = 1 contribution arises from a domain-scale 
meandering of the jet, termed "satellite modes" by Ref. 
[26] . The right panel of Fig. [5] shows that CE2 reproduces 
the contributions to the second cumulant at the longitu- 
dinal scales of the forcing, but is incapable of reproduc- 
ing the contribution from the satellite modes, when the 
system is this far from equilibrium. Interestingly these 
modes are also absent from quasilinear DNS calculations 
[6], which would seem to indicate that they arise as a 
result of eddy + eddy — >• eddy interactions. 

For systems driven even further out of equilibrium, 
CE2 struggles not only to reproduce all the structures 
of the second cumulant, but also the number of jets 
and their strength. As noted earlier, for smaller Rp the 
jets are more intermittent and meander more. Although 



zonal averages can be calculated, the constant meander- 
ing of the jets in latitude reduces the average jet strength. 
CE2 eventually settles down to a fixed point though we 
do not believe this to be a unique solution, with hystere- 
sis evident. The solution overestimates the strength of 
the jets and therefore the Rhines scale associated with 
them; hence CE2 has a tendency to underestimate the 
number of jets as shown in Figs. [3jb) and|4jc). 

This paper has demonstrated that DSS as approxi- 
mated by CE2 performs well in directly calculating the 
statistics for /3-plane turbulence in quasi-equilibrium. It 
confirms the earlier result [5j [6] that zonal jets do not 
require an inverse cascade to be driven, but can arise 
as the result of Reynolds stresses alone. However, and 
importantly, we have shown that as the system is re- 
moved further from equilibrium by reducing the zonos- 
trophy parameter i?^, CE2 can significantly overestimate 
jet strengths and predict the incorrect number of jets. 
We hypothesise that for such systems higher order cu- 
mulant expansions are required. If truncated at third 
order (CE3) the cumulant expansion includes eddy + 
eddy — > eddy interactions and should perform better in 
predicting statistics for out-of equilibrium systems. The 
potential utility of CE3 has been demonstrated for the 
problems of an isolated vortex [27] and fluid flow relax- 
ing to a prescribed jet [28]. We conclude by noting that 
although we have stressed the limitations of CE2, we be- 
lieve that the local /3-plane system driven stochastically 
is one of the stiffest tests of this method; it is very diffi- 
cult in both DNS and DSS to reach the quasi-equilibrium 
limit; although progress may be achieved utilising DSS 
implementing semi-implicit timestepping. Nevertheless 
CE2 provides a good qualitative description for systems 
where the important competing effects arise from inho- 
mogeneity, anisotropy, and turbulent fluctuations about 
a non-trivial basic state. Furthermore, CE2 should pro- 
vide a more accurate description for turbulence driven 
by instabilities; here the turbulence usually acts so as to 
'soften' the instability towards marginality driving the 
system back to quasi-equilibrium. In those cases CE2 
is able to capture the essential physical mechanisms of 
saturation. 
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